MODULE md1_utils

   USE md1_global_consts
   USE md1_global_types
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: md1_get_prim_batch

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE md1_get_prim_batch(basis,sh_pair,batch,np)

      IMPLICIT NONE
      TYPE(md1_basis),      INTENT(IN)  :: basis
      TYPE(md1_sh_pairs),   INTENT(IN)  :: sh_pair
      TYPE(md1_prim_batch), INTENT(OUT) :: batch(:)
      INTEGER(INTK),        INTENT(OUT) :: np

      REAL(REALK), PARAMETER :: thresh = -2.30258d+00 * 13.0d+00
      INTEGER(INTK) :: Ish, Jsh, I,J, Jhi
      REAL(REALK)   :: A(3),B(3), P(3), RAB(3)
      REAL(REALK)   :: ExpA, ExpB, ExpP, ExpPI, ExpAR2, R2AB, ExpKAB

      Ish = sh_pair%I
      Jsh = sh_pair%J
      A(:) = sh_pair%A(:)
      B(:) = sh_pair%B(:)
      RAB(:) = A(:) - B(:)
      R2AB = DOT_PRODUCT(RAB,RAB)

      np = 0
      DO I = basis%exp_id(Ish), basis%exp_id(Ish) + basis%nprim(Ish) - 1
         ExpA = basis%exp(I)
         ExpAR2 = ExpA * R2AB
         Jhi = basis%exp_id(Jsh) + basis%nprim(Jsh) - 1

         IF (Ish == Jsh) Jhi = I
         DO J = basis%exp_id(Jsh), Jhi
            ExpB = basis%exp(J)
            ExpP = ExpA + ExpB
            ExpPI = one / ExpP
            ExpKAB = - ExpAR2 * ExpB * ExpPI

            IF (ExpKAB >= thresh) THEN
               np = np + 1
               batch(np)%exp = ExpP
               batch(np)%KAB = EXP(ExpKAB)
               batch(np)%CAB = basis%C(I) * basis%C(J)
               IF ((Ish == Jsh) .AND. (I /= J)) THEN
                  batch(np)%CAB = two * batch(np)%CAB
               END IF
               P(:) = (ExpA*A(:) + ExpB*B(:)) * ExpPI
               batch(np)%P(:)  = P(:)
               batch(np)%PA(:) = P(:) - A(:)
               batch(np)%PB(:) = P(:) - B(:)
            END IF

         END DO
      END DO

   END SUBROUTINE md1_get_prim_batch

!-------------------------------------------------------------------------------

END MODULE md1_utils

